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ABSTRACT 

We study the physics of wave propagation in a weakly ionised plasma, as it applies 
to the formation of multifluid, MHD shock waves. We model the plasma as separate 
charged and neutral fluids which arc coupled by ion-neutral friction. At times much less 
than the ion-neutral drag time, the fluids are decoupled and so evolve independently. 
At later times, the evolution is determined by the large inertial mismatch between the 
charged and neutral particles. The neutral flow continues to evolve independently; the 
charged flow is driven by and slaved to the neutral flow by friction. We calculate this 
driven flow analytically by considering the special but realistic case where the charged 
fluid obeys linearized equations of motion. We carry out an extensive analysis of 
linear, driven, MHD waves. The physics of driven MHD waves is embodied in certain 
Green functions which describe wave propagation on short time scales, ambipolar 
diffusion on long time scales, and transitional behavior at intermediate times. By way 
of illustration, we give an approximate solution for the formation of a multifluid shock 
during the collision of two identical interstellar clouds. The collision produces forward- 
and reverse J shocks in the neutral fluid and a transient in the charged fluid. The latter 
rapidly evolves into a pair of magnetic precursors on the J shocks, wherein the ions 
undergo force free motion and the magnetic field grows monotonically with time. The 
flow appears to be self similar at the time when linear analysis ceases to be valid. 

Key words: diffusion — MHD — waves — shock waves — ISM: magnetic fields — 
ISM: clouds 



1 INTRODUCTION 

. It is well known that shock waves in weakly-ionised interstel- 
1 lar plasmas have a multifluid structure, where the charged 
and neutral compone nts of the p lasma behave as separate, 
interacting fluids ( MullanI Il97ll ). Because their m ultifluid 
natur e has profound observational consequences (|Drainel 
Il980l ). multifluid shocks have been the subject of numer- 
ous studies o n their structure, chem istry, and emission (last 
reviewed by iDraine fc McKee|[l993 ). Although the major- 
ity of these studies have assumed steady flow, a small frac- 
tion have carried out time dependent si mulations, either to 
follow the d evelopme nt of instabilities (iMac Low fc SmithI 
ll997l : IStone| [T997: No ufeld fc Stondll997l ') or to study evolu- 
tionary effects (|Smith fc Mac Low jl997|: C hiczc. Pineau des 
Forets fc Flower 1998: ICiolek fc Robergei 12002. : LesaflFre et 
al. 2004a,b). 

This paper describes the formation of a multifluid shock 
wave by a sudden disturbance, e.g., the collision of two 
cloud cores or the impact of a protostellar outflow onto sur- 
rounding core material. We focus on time scales < 100 rin, 
where the ion-neutral drag time, rin, is the slowing-down 



time for an ion drifting through a neutral gas. This is a very 
short time: rin ~ 0.01 yr in a typical dense core (^. Never- 
theless, there are good reasons for studying this extremely 
brief, unobservable period. First, the physics is interesting. 
The response of a plasma to a sudden disturbance depends 
on the wave modes it supports over a broad range in fre- 
quency, z/. In a weakly ionised plasma the allowed modes 
represent propagating waves if > and diffusion if 

V < Tj~^ ( ij3.3p . We wish to understand how these qualita- 
tively different behaviors manifest themselves at early times, 
and whether the resulting effects influence the flow at later, 
observable times. Second, the present study serves as a pro- 
totype for future work, on the effects of charged dust grains 
on m ultiffuid shocks. These effects are kno wn to be impor- 
tant l|Wardldll99ilCiolek fc Robergell2002|: Ciole k. Roberge 



& Mouschovias 2004; IChapman fc Wardk 



120061 ). Including 



dust will be analogous in some ways to the present study, 
but dust also adds new physics and a higher level of com- 
plexity to the problem. Third, the results of this paper have 
practical use. As a computational expedient, some time de- 
pendent simulations of multiflui d shock waves n eglect the 
inertia of the charged fluid (e.g., , Smith fc Mac L ow 1993). 
Since the inertia is important on precisely the time scales 
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studied here, the present work provides benchmark tests on 
this assumption. 

The plan of this paper is as follows. In ^we give the 
equations of motion for a time dependent, multifluid shock 
wave in a form which exploits the small time scale of interest. 
In 33]the linearized versions of these equations f ij3.1[l are pre- 
sented and we discuss the allowed wave modes (ii H3.2fl53)l . 
We emphasize that linear analysis yields highly accurate so- 
lutions for some special but realistic cases; the general case 
will be studied numerically in a separate paper. Analytical 
methods for calculating the time dependent flow of charged 
and neutral particles are described in ii H3.4ffX^ These tech- 
niques are used in |4]to find an approximate solution for the 
formation of a multifluid shock during the collision of two 
identical clouds. Our results are summarized in 9S] 



2 GOVERNING EQUATIONS 

We are interested in the multifluid flow which ensues when 
a weakly ionized plasma is accelerated and compressed by 
a sudden disturbance. We model the plasma as separate 
charged and neutral fluids which interact via elastic scatter- 
ing. The charged fluid is composed of ions and electrons (the 
effects of charged dust grains will be considered in a sepa- 
rate paper) plus a magnetic field which is everywhere frozen 
into the charged fluid. We adopt planar geometry with the 
magnetic field along the z direction and fiuid velocities along 
the ±x directions. 

The equations of motion for an arbitrary disturbance 
in a two-fiuid plasma were derived by Draine (1986; see 
also Nemirovsky, Fredkin & Ron 2002) . Here we solve modi- 
fied versions of these equations which take advantage of the 
extremely short time scale — typically less than 1 yr — on 
which we follow the flow. The charged fluid is described by 
the equations of mass conservation, momentum conserva- 
tion, and the induction equation: 



dpi d , , 



0, 



'dt 

and 

dB 



+ Vi- 



dvi 



Pi dx 







(1) 



(2) 



(3) 



respectively, where pi and Vi are the density and velocity of 
the charged fluid, Vn is the velocity of the neutral fluid, and 
B is the magnetic field. The first term on the RHS of eq. ^ 
is the frictional acceleration produced by elastic scattering 
between ions and neutral particles and the second is the 
acceleration caused by the magnetic pressure gradient. 

The characteristic time scale for acceleration by friction 
is the ion-neutral drag time, rin. If the charged and neutral 
fluids were each composed of a single species, then 



1 + mi / m„ 



(4) 



where rrin and mi are the neutral and ion particle masses, 
respectively, Un is the number density of the neutral fluid, 
and is the momentum transfer rate coefficient for 

elastic ion-neutral scattering. For a typical cloud core with 



Un = 2 X 10'* cm~^ and mi = 25mn, one finds rin ^ 0.01 yr. 
This is a fundamental time scale for the flow. 

We have omitted the energy equation for the charged 
fluid because eq. (HI-lIS} are independent of the ion and elec- 
tron temperatures, Ti and Tc. This is appropriate at the 
very low fractional ionizations 10~*) of interest here, 
where the pressure of the charged fluid is dominated by 
magnetic pressure. The cooling rate of the neutral fluid gen- 
erally depends on Ti and Tc (via the rates of ion- and elec- 
tron impact processes) but we assume that radiative cooling 
is negligible on the time scales of interest. For a gas with 
Wn = 2 X 10''cm~^, the cooling time is > lyr if the neu- 
tral temperature is less than ~ 2300 K. In eq. ((l| we have 
omitted a term which represents mass transfer between the 
charged and neutral fluids. This is always a good approxima- 
tion because the recombination time scale is always ^ 1 yr 
in dense clouds. 

We assume that the neutral fluid is governed by Euler's 
equations for adiabatic flow, 

dpn d . , 
_ + _(p„„„)^0. 



dVn _|_ ^ dVn 

dt " dx 
and 

Pn=KpZ. 



1 dPn 
pn dx ' 



(5) 
(6) 

(7) 



Mass transfer in eq. ([Sj is also neglected, for reasons noted 
above. We have neglected momentum transfer by friction in 
eq. (|6]) because the time for friction to accelerate the neutral 
fluid is 

Tni = ^ Tin, (8) 
Pi 

or about 10* yr in a dense core0 We have also neglected 
heating of the neutral fluid by ion-neutral friction and the 
associated acceleration by thermal pressure gradients; one 
can show that these effects are small on time scales ~ 1 yr. 

The equations of motion for the charged fluid depend on 
the neutral velocity, Vn, and number density, Un- In the next 
section, on the dynamics of the charged fluid, we assume 
that Vn and rin are known functions of x and t which have 
been determined by solving Euler's equations. How this all 
works out for a particular example is demonstrated in ^ 



3 THE PHYSICS OF DRIVEN WAVES 

3.1 Linearized Equations of Motion 

We follow the flow of the charged fluid by solving the lin- 
earized versions of eq. ©-(IS}. This is an expedient which 
allows us to work the problem analytically. Our analyti- 
cal solutions reveal the essential physics; numerical solu- 
tions of the full nonlinear equations will be discussed else- 
where. However it is important to note that linear theory 
is highly accurate for some initial conditions at sufficiently 
short times; this is ultimately because the speed of a typ- 
ical disturbance (~ 10-100 km s~^) is much smaller than 



^ Of course this means that the total momentum of the charged 
plus neutral fluids is not conserved but the associated error is 

0{nn/rni) ~ 10-6. 
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the ion Alfven speed (~ 1000 km s^^). A specific example is 
discussed in 311 

The zero-order solution of eq. ([l|-(|3]) is a spatially uni- 
form state in which the charged fluid has constant density 
pio, magnetic field, Bq, and zero velocity. Onto the zero- 
order state we superpose density, magnetic field, and veloc- 
ity perturbations denoted r, 6, and u, respectively. We adopt 
dimensionless perturbations so that 



piix,t) 
B(x,t) 
and 



Pio [l + r{x,t)], 
Bo [l + b{x,t)], 



Vi{x,t) = ViAO u{x,t), 



where 



Bo 



ViAO 



\/47rpio 



(9) 
(10) 

(11) 
(12) 



is the ion Alfven speed in the zero-order state. We also adopt 
dimensionless time and distance units. In all subsequent dis- 
cussion, t is dimensionless time in units of Tin and x is di- 
mensionless distance in units of iiiAoTin- In this paper we 
shall assume that Tin is independent of x and t. This allows 
us to obtain analytical solutions but, because Tin depends 
on the density of the neutral fiuid, only hypothetical fiows 
with uniform density can be studied analytically. 

Linearizing eqs. (HJ-Q about the zero-order solution 
gives the equations of motion: 



r + u' = 0, 
b + u' = 0, 
and 

il + U + b' = Un, 

where 

Un {X,t) = Vnix,t)/ViAO 



(13) 
(14) 

(15) 

(16) 



and the dots and primes denote partial derivatives with re- 
spect to t and X, respectively. We seek solutions of equations 
(fT3)) - (fT5l) subject to the initial conditions 

(17) 
(18) 



r(x,0) = ro{x), 
b{x,0) = bo{x), 
and 

u{x, 0) = 1*0(2;), 



(19) 



where ro, bo, and wo are small but otherwise arbitrary func- 
tions. We assume that Un has been determined as described 
in 3U For the purposes of this section it is a known, small, 
but otherwise arbitrary function. 



3.2 Fourier Analysis 

To proceed we Fourier transform eq. (|13|l - (|15p to eliminate 
the spatial derivatives. We define the Fourier transform of 
fix,t) by 



fik,t): 



1 

2^ 



Writing the perturbations as Fourier integrals transforms 
eq. (|13|) - (|15|) into a set of linear, inhomogeneous, coupled, 
ODEs for the transforms of the perturbations: 

dy 



dt 



= Qy + f, 



where the vector of unknowns is 
y{k,t) = [f{k,t),h{k,t),u{k,t)Y , 
the couphng matrix is 



Q(fc) 



and the source term, 

i{k,t) = [0,0,Un(fc,t)]*, 




(21) 



(22) 



(23) 



(24) 



represents the effects of frictional driving. 

Equation (|2ip can be solved by standard methods. The 
solution has the form 

y(fc,t) =y„(fc,t) + yp(fc,t), (25) 

where the particular solution, yp, is any solution of eq. (|21|l 
and the homogeneous solution, yj^, is any solution with f = 
0. The "total" solution must satisfy the initial conditions 

yp(fc, 0) + y,,(fc, 0) = [f 0, bo,uo\ * (26) 

but there is a degree of arbitrariness in how the right side of 
(|26|) is apportioned between the particular and homogeneous 
parts on the left. We require 

yp(fc,0) = (27) 
and 



yhC^.O) = [fo,bo,uo\^ 



(28) 



Then for a given set of initial conditions, the homogeneous 
solution represents the disturbance that would occur, for 
the same initial conditions, if the driving force was zero. We 
show in >j3.3l that the homogeneous solution has no growing 
modes; it is therefore entirely transient in nature. The par- 
ticular solution always has a transient component but may 
also contain a steady part if the driving is steady (see tj3.9p . 

3.3 Dispersion Relations 

The solution of eq. (|2ip depends on the eigenvalues and 
eigenvectors of the coupling matrix. It is easy to show that 
Q has three eigenvalues, — --luJOi and —luj+t where 



LOO = 0, 

and 

oj±{k) = -'-±R{k). 

The corresponding eigenvectors are 

lo = [1,0,0]* 

and 

i± = [k/uj±,k/u)±,l]^ , 



(29) 
(30) 

(31) 
(32) 



dx e 



■'fix,t). 



(20) 



respectively. 

The eigenvalues and eigenvectors depend on the 
complex-valued function 
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Irak 



Rek 



Figure 1. Function R{k) has a branch cut in the complex k plane 
along the interval indicated by the sawtooth curve. 



R{k) = {k^ ~ k^,) ^ , (33) 

where kc — 1/2 is the dimensionles^ "critical wave num- 
ber." For real k the meaning of R is ambiguous on the in- 
terval — fee < k < +kc, where R has a branch cut in the 
complex k plane (Fig. [1} . We will take 



R{k) 



k < -kc 



k > +kc 



(34) 



where ".y " denotes the positive square root. This means 
that for real k we evaluate R{k) "just above" the branch 
cut, a fact which becomes crucial when Fourier transforms 
are evaluated as contour integrals (App.|A])- With our defi- 
nition of R, the dependence of ui- and lo+ on k (for real k) 
is as shown in Fig. [2]^ Kulsrud and Pearce (1969) found 
the analogous dispersion relations for Alfven waves in a 
cold io n-neut ral plasma. For dimensionless wave numbers 
k ^ •\/ Pi/pn ~ 10""^, where the neutrals act as a stationary 
background, th eir dispersion relation re duces to expression 
pOl) for Lo± [cf. iKulsrud fc Pearce|[l969l . eq. (C7)]. 

The physics of the different wave modes is well under- 
stood but worth repeating for later discussion. The loq mode 
represents perturbations (e.g., perturbations in B / pi) that 
move along with the charged fluid. For |fc| > fcc the lj+ and 
u- modes represent propagating waves (Fig. For very 
large wave numbers, the dispersion relations become 



lim uo± — ± fc. 

|fc|»i 2 



(35) 



In this limit the phase velocity approaches ±1 (^ ±UiAo 
in ordinary units) and the damping time approaches 2 (^ 

For |fc| < fee the ui+ and ui- modes are both evanes- 
cent (Fig. [2]). However they describe qualitatively different 
behavior. For very small wave numbers the dispersion rela- 
tions become 

^ In ordinary units, kc = 2n/\c, where Ac = 47r»;iAof"in is 
the maximum wavelength for propagating ion magnetosound 



(and Alfven) wave s in a cold plasma (Kulsrud . 
ICiolek. et al1l2004l') . 



Pcarcg [l969t 



Pi 




-1.0 



Figure 2. Real part of tj_|_ (solid) and ui^ (dashed) plotted vs. 
wave number k. The curves are degenerate for |fc| ^ 1/2. All 
quantities are dimensionless. 
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Figure 3. Imaginary part of aj+ (solid) and ui— (dashed) plotted 
vs. wave number k. The curves are degenerate for |fc| 1/2. All 
quantities are dimensionless. 



Q)_ 

I I I 



lim uj+{k) — —ik . 
and 

lim (k) — —1 + ik^ 



(36) 



(37) 



These are to be compared with the dispersion relation 
<^diff(fc) = —ik^oi (38) 

for the diffusion equation with diffusion coefficient a. Evi- 
dently u+ represents diffusion with oi = 1 t'?Ao'''in) and 
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zero damping. This "difFusion mode" dominates all solutions 
at large times. In contrast, aj_ represents "antidiffusion" 
with a = —1 and strong damping. The u- mode describes 
the transient compression of B which can occur at early 
times for some initial conditions. However it is never impor- 
tant for t > 1. 



3.4 Homogeneous Solution 

The homogeneous solution of eq. (|21|l is 



(39) 

where summation is implied by the repeated index. The co- 
efficients Am are functions of k but not time and are deter- 
mined by the initial conditions. Setting t = in eq. (I39|l and 
substituting the result into the LHS of eq. (|28p gives a set 
of linear algebraic equations for {Am}. The solution is 



A 



2R\ k 

Ao = fo~- bo, 
and 



•0 — "U-O 



A 



+ 



2R\ k 



(40) 
(41) 

(42) 



The expansion coefficients obey the symmetry requirement 

A- —> A+ under the interchange > To see this, note 

that R={iu+- L0-)/2 [cf. eq. pO]! ]. 

The transforms of the perturbations are obtained by 
substituting the expansion coefficients into eq. (|39p . After 
lengthy algebra we find that 



fh(fe, i) = 6h + ^0 - bo, 

feh (fc, i) = 2ttG bo + 2Tvd bo - 2ttG' iio , 
and 

jih {k,t) = 2tvG uo - 2ttG' bo . 
We have introduced the function 
G{k,t) = G+{k,t) - G-{k,t), 
where 

G^{k,t)^^' 



Atv R 



and 



G-{k,t) 



47r R 



(43) 
(44) 

(45) 

(46) 

(47) 
(48) 



G'(k, t) and G{k, t) are the Fourier transforms of G' {x, t) and 

G{x,t), respectively. We shall see shortly that G, G' , and G 
are the Fourier transforms of certain Green functions. The 
Green functions are calculated in Appendix]^ and discussed 
in i )3.6l Here we simply note that they are causal: 

(49) 



G{x, t) = G{x, t) = G'(x, t)=0 if |x| > t. 



The homogeneous solution for the perturbations is ob- 
tained by taking the inverse Fourier transforms of expres- 
sions (|43|l - (|45p . This can be done by inspection using the 
Convolution Theorem, 



dke^"- fik)g{k,t) 



1 

2^ 



The result is the homogeneous solution: 
rh(x,t) = bi,(x,t) +ro{x) - bo{x) 
6h(x, t) = {G\bo) + {G\bo)- {G' \uo ) 
and 

ui,{x, t) ^ (^G\uo) - {G' \bo ) , 

where the angle brackets signify convolution. 



dx' g {x — x' ,t) f {x') , 



(51) 
(52) 

(53) 
(54) 



and causality has been used to refine the limits of integra- 
tion. 

Expression (|5ip for the density perturbation has a sim- 
ple physical interpretation. Taking the ratio of eq. (|10|) and 
© and neglecting terms of second order gives 



BiX, t) Bo , , , , . , 

^ ' ^ - — [l + &(x,t)-r(x,t)]. 



(55) 



pi(x,t) pio 

But according to eq. (|51|) . the factor in square brackets is a 
conserved quantity so 

B{x,t) _ B{x,0) 



pi{x,t) pi{x,Q)' 



(56) 



Noting that = in the unperturbed state, we see that 
expression (|56p implies fiux freezing. Of course flux freezing 
was put into the solution a priori, so there was really no 
need to calculate the density perturbation independently of 
the magnetic field perturbation. However it is reassuring to 
know that the homogeneous solution is self consistent in this 
respect. 



3.5 Particular Solution 

The particular solution of eq. (|21|l has the form 



yp(fc,t) = Gme 



(57) 



where the coefficients Cm generally depend on time as well 
as k. One can verify that expression (|57|) is a solution of 
eq. (I21|l provided 



Cme '""'{^ = t{k,t). 



(58) 



Expression (|58p is a set of coupled ODEs for {Cm} with 
initial conditions 



Cm(fc,0) =0. 
The solution is 

C-(fc,t) = / dt' Un{k,t') 



Co{k,t)=0, 
and 



C+ik,t) 



2R 



dt' Unik,t') e"+*' 



(59) 

(60) 
(61) 

(62) 



dx' fix') g{x-x',t).{50) 



Substituting the expansion coefficients into expression 
H57|) yields the Fourier transforms of the perturbations. The 
inverse transform then yields the particular solution: 



rp{x,t) = bp{x,t), 



(63) 



6 W. G. Roberge and Glenn E. Giolek 



35 I ■ III I ■ III I ■ III I ■ III I ■ III I III ■ 



.30 : 
.25 : 
C^ .20- 
O .15^ 



.10 - 



.05 - 



.00 



A 




t=l 



-15 -10 -5 5 10 15 



Q5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



X 

< 



.00 



-.05 



-.10 



-.15 




t= 1 



I ■ ■ ■ ■ I ■ ■ ■ ■ I ■ I ■ ■ I ■ ■ ■ ■ I 



-15 -10 -5 5 10 15 



Figure 4. Green function G{x,t) plotted vs. x at t = 1, 5, and 
10. Note the discontinuities at x = itt. 



bp{x,t) 



and 



■ip(x,t) = 



dt' I dx' Un (a;', i') 

J x-t + t' 

G' [x x' ,t ~ t'^ , 



t fX+t—t 

dt' / dx' Un (x', t') 

J x-t + t' 



X G{x-x',t-t') , 



(64) 



(65) 



where the Convolution Theorem has been used again. Com- 
paring eq. (|63p to eq. (ISlfl and noting that the particular so- 
lution vanishes at t = 0, we see that the particular solution 
also incorporates flux freezing. Since flux freezing uniquely 
relates the density perturbation to the magnetic field per- 
turbation, we omit further discussion of r{x,t). 



Figure 5. Function A{x,t) plotted vs. x for t = 1, 5, and 10. 
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3.6 Green Functions 

The Green function G is calculated in Appendix |X] We find 
that 

G(x,f) = ie-'/^/o(C/2)n(j), (66) 
where Iq it the modified Bessel function of order zero, 
i = ^f^-x^, (67) 
and the rectangle function, 

r 1 if l^l < 1/2 
= \ , (68) 

[ if \x\ > 1/2 

insures causality. In Fig.|4]we plot G vs. x at selected times. 
Differentiating G yields the other Green functions, 



Figure 6. Function r{x,t) plotted vs. x for t = 1, 5, and 10. 

G{x, t) = A(x, t) + i e-'^^Six + t) + ^ e~'^^S{x - t) (69) 
and 

G'{x, t) = T{x, t) + i e-'''^S{x -f- 1) - i e-*''^5{x - t), (70) 
where 



k(x,t) = ~^G{x,t) + ^ie"'''^ 



h (e/2) 



n(£) (^1) 



and 



V(x,t) = --^xe-'/^ 



h (C/2) 



(72) 
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and 



Plots of A and T are given in Fig. [SHS] 

It is useful to know the normalization of the Green func- 
tions (e.g., to check the numerical evaluation of convolution 
integrals). It is easy to show that 

dx G{x,t) = 1- e~\ (73) 
dxG{x,t)=e~\ (74) 

dxG'{x,t)^0. (75) 
3.7 Large Time Behavior 

The discussion of the dispersion relations ( H3.3|l suggests 
that diffusion dominates the physics on large length scales. 
Because the damping rate of the diffusion mode increases 
with k (Fig. [31), this should be reflected in the large-time 
behavior of the Green functions. In Appendix [B] we show 
that 

1 



G{x,t) — > Gdia{x,t) 



/47rt 



as t 



(76) 



Consistent with the discussion in i )3.3l G{x,t) approaches 
the Green function for the diffusion equation with unit dif- 
fusion coefficient. We also show in App. |B] that the other 
Green functions are smaller than G at large times, with 

-1/2 



G:G':G 



1 : r 



t 



(77) 



(see Fig. [T]). 

The asymptotic forms of the Green functions lead to a 
particularly simple result for the solution when driving is 
absent. If one replaces G with Gditr and uses expression (|77p 
to omit all but the nonvanishing terms of largest order in 
the homogeneous solution, the latter becomes 



b^{x,t) ^ (Gdift l&o) 
and 

UY,[x,t) ^ - (Gdifl- \bo) . 

Equation (|78|l indicates that the magnetic field satisfies a 
diffusion equation when t ^ 1, as expected. Equation (|79[) 
is understood by noting that 



(78) 
(79) 



d 



(Gdiff- 1^0 ) = ^ (Gdift' \bo } 



so that eq. (|79|l is the same thing as 
Uh{x,t) ^ —b'i-,{x,t). 



(80) 



(81) 



The last expression says that the magnetic and drag forces 
on the charged ffuid balance one another [cf. eq. p5|l with 
Un = 0]. This is also expected: when f S> 1, the inertia of the 
ions becomes negligible and they undergo force-free motion. 



3.8 Small Time Behavior 

When t <^ 1, the effects of friction are small and the results 
of ^ for driven waves should reduce to the corresponding 
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Figure 7. The functions G, t^/^G', and tG plotted vs. X for 
t = 100. 



results for ideal MHD. The latter can be found in a form 
suitable for comparison by setting m = itn = in the lin- 
earized momentum equation and retracing the steps leading 
to eq. (f5T|| - (f53)) (for the homogeneous solution) and eq. ((63])- 
(|65l) (particular solution). The calculation is very straight- 
forward and we simply give the results: 

bnf{x,t) = (Gnf \bo) - (Gnf | UQ ) , 

and 

Unl{x,t) = (Gnf \uo ) - (Gnf 1 60 ) , 

where the "frictionless Green functions" are 

Gnf(x,t) = i 5(X -f t) -f i 6ix - t) 

and 

Gnf(x,t) = \5{x + t)~\ 5{x ~ t). 



(82) 



(83) 



(84) 



(85) 

We wish to compare expressions (|82|) and (|83|) to the 
corresponding solution with friction in the limit t <^1. The 
particular solution obviously goes to zero as f — > 0. To find 
the limit of the homogeneous solution, we note that any 
convolution integral approaches zero if its kernel is G, A, or 
r, because the amplitude of each of these functions remains 
finite as t ^ Q, and is zero for |x| > t. Retaining only the 
delta function kernels (with exp {—t/2) ~ 1 for t <^ 1) is 
equivalent to making the replacements 



(86) 
(87) 



G{x,t) 0, 
G{x,t) Gni{x,t), 
and 

G'(x,t)^GUx,t). 

in eq. (|52p - (|53|) . This shows that the physics of driven waves 
reduces to ideal MHD in the appropriate limit. 

A comparison of the Green functions for waves with and 
without friction also shows how ion-neutral friction alters 
the physics of wave propagation. This is best illustrated by 
examining the solutions in terms of the Riemann invariants, 
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Figure 8. In ideal MHD the Riemann invariants J± are con- 
served quantities which are transported along characteristics. In 
the presence of friction they are transported along the same char- 
acteristics but not conserved. 



J_ {x, t) = b{x, t) — u{x, t) 
and 



(89) 



J+{x,t) = b{x,t) +u{x,t). (90) 
In terms of J_ and J+ the solution for ideal MHD becomes 



J-{x,t) = J_ (x + t,0) 
and 

J+{x,t) = J+{x- t,0) , 

where we have used the fact that 

{f\5{x±t)) = f{x±t) 



(91) 



(92) 



(93) 



to evaluate the convolutions. The analogous solution for 
driven waves is 

J-{x,t) = J_ {x + t,0) e''^''+{G\bo)+{A\J-^)+{r\J+) (M) 
and 

J+{x,t) = J+ {x-t,0) e-*'''+(G|6o)-(r|J+)+(A|J_).(95) 

Expressions (|9ip and (|92p say that the Riemann invari- 
ants are conserved quantities in ideal MHD. The former 
is transported along the characteristic curve dx/dt = —1 
and the latter along dx/dt = -1-1 (Fig. [S]!. It follows that 
the solution at each spacetime point (x, t) depends only on 
the initial conditions at the points {x — t, 0) and {x + t, 0) 
(Fig. [8]). This is the essence of wave propagation. When fric- 
tion is present, the situation is similar yet not identical. The 
first term on the RHS of eq. (I94|l and eq. (I95|l may be in- 
terpreted to mean that the Riemann "invariants" are still 
transported along characteristics, but with exponential at- 
tenuation. That is, they represent wave propagation with 
damping. However the presence of additional terms shows 
that this is not the whole story. For example, the nonzero 
widths of the convolution kernels on each RHS imply that 
the solution at (a;, t) depends on the initial conditions over 
the entire spatial interval {x — t, x + t). This is the signature 
of diffusion. 

The physics of wave propagation, diffusion, and the 
transition from wavelike to diffusive behavior is apparent in 



the plots of the Green functions (Fig.|4H6]). Wave fronts are 
located by the jumps at a; = ±t, propagating with velocities 
of ±1. Diffusion is implied by the finite widths of the Green 
functions. The transition from wave propagation to diffusion 
occurs at times of order unity, when the jumps have decayed 
conspicuously. At much greater times the jumps become in- 
visible and all vestiges of wave propagation disappear. 



3.9 Steady Driving 

The physics of driven waves simplifies dramatically when the 
frictional driving is steady. Since we are interested in time 
scales where the neutral flow is approximately steady, it is 
worth discussing this case in some detail. Only the particular 
solution needs to be considered; the homogeneous solution 
depends only on the initial conditions, and so requires no 
modification. 

If the neutral flow is steady then it is possible to write 



-i(x, t) — D {x — u^t) . 



(96) 



where D(x) is the neutral velocity profile and itg is the pat- 
tern velocity0 Then the Fourier transform of Un has the 
form 



Un{k,t) = D(k)e- 



(97) 



and the simple time dependence of u-n makes it possible to 
find the expansion coefficients Cm explicitly. Substituting 
expression (|97|l into equations (|60p and (|62p and evaluating 
the integrals gives 



C-e"'"'-' =b{k)S-{k) [e 
and 

C+e-'"+* = £)(fc)S'+(fc) [e 
where 



Q(k) 



2k — lUf, 



A: (1 + m|) — jUg 



2k — iu„ 



Pik) = 2 



k{l + u'g) — lUg 



(98) 
(99) 

(100) 
(101) 

(102) 



Having found C_ and C+, we only need to calculate 
the inverse transform of expression (I57p to find the particu- 
lar solution. Evidently the latter contains two parts: terms 
in C± that are proportional to exp {—ikugt) represent steady 
fiow of the charged fiuid, which must eventually result from 
steady driving by the neutrals. The other terms in C± rep- 
resent the transient response to driving. 

First we evaluate the velocity, Up. Substituting expres- 
sions ((98)) and ((99)) into eq. (|57)) gives 

^ip = + Upt, (103) 
where 



Note that the linearized equations of motion arc not invariant 
under Galilean transformations: they are valid only in the frame 
where the undisturbed gas far upstream is at rest. The pattern 
velocity is evaluated in this frame. 
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«p, = D {S-+ S+) e-'''^ 
and 



Upt 



(104) 



(105) 



are the steady and transient parts, respectively. Now the 
inverse transform of expression (|104|l is 



where 

S{k) = S-+S+ = 



k — tu„ 



(106) 



(107) 



The inverse transform of S is evaluated by simple contour 
integration to find 

{0 if u^x < 

(108) 
2tt \ug\ exp {—Ugx) if UgX > 

If we use this result and the Convolution Theorem to eval- 
uate expression p06|) . we obtain the steady part of the ion 
velocity: 

Ups{x,t) = D{x-Ugt), (109) 
where 

/' + 00 



D(x) 



2n 



dx' S (x') D (a; — x'^ 



(110) 



This is a remarkably simple result. If the neutral flow is 
steady, then the charged fluid will eventually undergo steady 
flow with the same pattern velocity, consistent with com- 
mon sense. When both flows are steady, the velocity of the 
charged fluid is just the convolution of the driving force, D, 
with an exponential response function, S. 

Now consider the transient part of the ion velocity. Sub- 
stituting expressions (|100|) - (|102p into equation ([105}, we 
find after routine algebra that 



Upt{k,t) = 2ty 



[2iQ -P)DG- 2PDG 



(111) 



Inverting the transform is performed by a twofold applica- 
tion of the Convolution Theorem once the inverse transforms 
of Q and P are known. We find that 



P{x) = ]^S{x) 
and 

Q{x) 



■S'{x) - ^S{x). 



(112) 



(113) 



2?tg ^ ' 4 

Using these results and the Convolution Theorem yields 
upt{x,t) = {G\D) - {g\d) - {g\d). (114) 
The transform of the magnetic field perturbation is 

bp^—C-e-"^-' + —C+e-"^+\ (115) 

The inverse transforms are obtained by steps very similar 
to those used to find the velocity. We omit the details and 
simply state the result: 

bp{x,t) = D {x — Ugt) /ug — (g\D /ug^ 

-{G\b/ug) + {G\Dug). (116) 



4 A CLOUD-CLOUD COLLISION 

As an application of the methods developed in ^ we con- 
sider the collision of two identical, uniform, semi-infinite 
clouds. The relative motion is taken to be along x and the 
magnetic fields to be everywhere along z. Prior to the col- 
lision the free surface of each cloud is a plane normal to 
X. The collision occurs at t = 0, when the free surfaces 
touch at the contact discontinuity, a; = 0. We take the 
relative velocity of the clouds to be An = 20 km s~^ and 
adopt initial conditions appropriate for a dense core: the 
neutral fluid in each cloud is pure molecular hydrogen with 
Tin = 2 X lO^cm"^, the charged fluid has mi = 25amu, 
rii/n-n = 3 X 10^®, and the unperturbed magnetic field 
is Bo = 50 iiG. Then the characteristic speeds, time- and 
length scales are v\ao = 894 km s"'^, rin = 1.1 x 10~^ yr, 
Tni = 3.0 x lO"* yr, and wiao rin = 3.2 x 10^^ cm. 

The homogeneous solution depends only on the initial 
conditions, which are 



ro(x) = bo{x) = 



(117) 



for the density and field perturbations. We adopt a refer- 
ence frame (the "CM frame") where the clouds approach 
with equal and opposite velocities. In this frame the initial 
velocity perturbation is 



uo{x) = -iAusgn(a;), 
where 

Au = — .022 



ViAO 



and 



sgn(a;) 



-1 x<Q 



+1 x>Q 



(118) 



(119) 



(120) 



Using these initial conditions, we calculated the homoge- 
neous solution by evaluating expressions (|52p and (|53|l FI 

The particular solution depends on the velocity profile 
of the neutral fiuid, which is determined by solving Euler's 
equations. This is straightforward as the example considered 
here is just a special case of the Rie mann proble m in gas 
dynamics. Using the Riemann solver of iTorol l] 19991 ). we cal- 
culated Vn{x,i) for an adiabatic ideal gas with 7 = 5/3. We 
set the temperature of the unperturbed neutral gas some- 
what arbitrarily to Tn^o = 10 K; that is, we did not find 
Tn.o self-consistently by requiring thermal balance in the un- 
perturbed state. (Since the shocks in this example are very 
strong, the value of T'n,o hardly affects the solution.) The 
neutral flow is steady. It consists of forward- and reverse J 
shocks propagating with constant velocities in the ±x direc- 
tions. The shock velocities are ±3.3 km s~^ in the CM frame, 
corresponding to identical shock speeds of 13.3 km s~^ rela- 
tive to the upstream gas. Having found v-n{x, t), we obtained 
the particular solution from expressions (|109|) . (|114|l . and 
(fTT6)) . 

The velocity of the neutral fluid is plotted in Fig. [5] at 
the extremely early time t — 0.001 (^ 300s!). The region 
between the J fronts (the "g-star region") contains neutral 



Taking care to evaluate expressions I I52I I and 11531 in reference 
frames where they are valid. See footnote [3] 
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Figure 9. Velocity of the neutral fluid in the CM frame at the 
extremely early time t = 0.001. The discontinuities are forward 
{x > 0) and reverse (x < 0) J shocks. They propagate away 
from the contact discontinuity, a; = 0, at a constant speed of 
3.3km s~^. 
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Figure 10. Velocity profiles of the charged (solid) and neutral 
(dashed) fluids at the extremely early time t = 0.001. The g-star 
region is only fs; 10* cm wide at this time and so unresolved. 
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Figure 11. Magnetic field perturbation AB = B — Bo at the 
extremely early time t = 0.001. 



gas which has been swept up, heated, and compressed by 
the shocks. Because these shocks are very strong (with Mach 
numbers > 50), the swept-up gas is denser than the undis- 
turbed gas by a large factor (« 4), in gross violation of our 
assumption that rin is constant. Moreover the g-star region 
is so hot (Tn « 8, 000 K) that radiative coohng is important 
even on time scales < 1 yr. In a separate paper we describe 
numerical calculations that allow for variations in rin, ra- 
diative cooling, and other (e.g., nonlinear) effects. Here we 
temporarily forego these complications and simply warn the 
reader that our example may not be correct in detail. 

The velocity profile of the charged fluid at t — 0.001 
is plotted in Fig. 1101 At this time the effects of friction are 
negligible, so the flow in Fig. [TO] is just the magnetohydro- 
dynamic analog of the flow in Fig. [9] The MHD flow has a 
growing "m-star region" bounded on either side by a moving 
J front. In both the gas dynamic and MHD flows, a particle 
of fluid passing through a J front is brought to rest by an 
impulsive force inside the front; the force is collisional in the 
first case and magnetic in the second. In the gas dynamic 
fiow, the collisional force vanishes outside the J fronts so the 
neutral fluid remains at rest inside the g-star region (Fig.|9]). 
In the MHD flow, the magnetic force vanishes outside the J 
fronts but friction does not. The charged fluid appears to be 
at rest in the m-star region (Fig. llOp only because t <^ Tin. 

The magnetic field at t = 0.001 is plotted in Fig. [11] 
The m-star region is conspicuous as a "magnetic core" of 
compressed fiuid and magnetic field centered on the con- 
tact discontinuity. The compression is approximately uni- 
form because the velocity gradient in the m-star region is 
approximately zero. The compression ratio is small {~ 1.01) 
because the "shocks" in the charged fluid have ion Alfven 
Mach numbers of only ~ 0.01. Consistent with the linear 
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Figure 12. Velocity profiles of the charged (solid) and neutral 
(dashed) fluids at t = 0.1. The horizontal scale is too coarse to 
resolve the g-star region. 
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Figure 13. As in Fig. I12l but plotted on a scale that resolves the 
g-star region. 



analysis, the velocity jumps {Svi = ±10 km s ^) and mag- 
netic field jumps {SB = 0.562 fiG) obey the relation 

M = M (121) 

Bo ViAO ^ ' 

for linear magnetosound waves. We note also that the math- 
ematical approach we have adopted treats discontinuities ex- 
actly (i.e., without smoothing) and that discontinuties prop- 
agate at the correct speeds. 



PQ 
<1 



.7 
.6 
.5 
.4 
.3 
.2 
.1 



I I I 11 I I 1 I I 1 I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

t = 0.1 



■ ■ 



.0 

-4 -2 











X (10^^ cm) 



Figure 14. Solid: Magnetic field perturbation at t = 0.1. The 
horizontal scale is too coarse to resolve the g-star region. Dashed: 
Result of an identical calculation with friction turned off. Notice 
that the area under each curve is the same. 



Figures [T2Ml4l describe the flow at t = 0.1 («-+ 3 x lO" s), 
when the J fronts are at j: = ±1.07 x 10^" cm in the neutral 
fluid and x = ±2.86 x 10^ cm in the charged fluid. Depar- 
tures from ideal MHD are now visible as slight reductions 
in 5vi (to ^ ±9.5kms-i) and 5B (to « 0.536 AtG). The 
ions inside the m-star region are moving. They flow toward 
the contact discontinuity with vi ~ ±0.5 km s^^ just down- 
stream from the J fronts, decreasing to i;i = at x = 
(Fig. [13}. The compression inside the m-star region is no 
longer uniform; however the density /magnetic field gradient 
is almost constant, with B and p increasing toward x = Q 
(Fig. [13. 

The qualitative properties of the fiow at t — 0.1 follow 
from simple physics. To see this, it is useful to note that a 
particle of charged fluid moving with Vi — 0.5 km s^^ would 
travel only ~ 10^ cm in 3 x 10* s. Each point on the solid 
curves in Fig. I12H14I therefore labels a fluid particle which 
has not moved appreciably since t — 0. Now consider the 
histories of various fluid particles, starting from the initial 
state t = 0.001 when the effects of friction on the flow were 
negligible. More specifically, consider the forces at t = 0.001 
on particles that are inside the m-star region excluding the 
J fronts. All of these particles have Wi = at i = 0.001 
(Fig. llOfl and the magnetic force on each vanishes (Fig. II If) . 
The friction force is small inside the g-star region {vn — Ui ~ 
0) but not outside (wn — Vi « ±10 km s"'^). Since the net 
force is small inside the g-star region, we expect the charged 
fiuid there to be almost at rest at later times. This is true at 
t = 0.1 (Fig. I13|l . The net (=frictional) force is much larger 
outside the g-star region; the sign of Vn — Vi is such that 
the particles tend to speed up and flow toward the contact 
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Figure 15. Magnetic, drag, and net forces in dimcnsionless units 
at t=0.1. Only points between the forward shocks in the neutral 
and charged fluids are shown. 
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Figure 16. As in Fig. 1151 but on a scale that resolves the g-star 
region. Note force scale. 



discontinuity. This explains the nonzero values and sign of 
Wi at f = 0.1 (Fig.[T2j. 

However the preceding argument does not explain the 
velocity gradient in Fig. 1121 which would have the opposite 
sign if friction was the only force at work. This is because 
particles close to a:: = have been accelerating longer than 



particles that have just emerged from the J front. The argu- 
ment is incomplete because, contrary to intuition, the mag- 
netic force is comparable to friction ai t — 0.1. Indeed, the 
two forces differ by less than ~ 1% (Fig. [15}. The magnetic 
force is caused by ambipolar diffusion. The acceleration pro- 
duced initially by friction produces ion motions that trans- 
port magnetic field lines toward the contact discontinuity. 
Now the total magnetic fiux threading the m-star region is 
the same, at any given time, whether friction is present or 
not (Fig. I14|) FI The motion of field lines toward x = must 
therefore leave a deficit of field lines farther out (relative to 
the frictionless case; compare the solid and dashed curves in 
Fig. I14|) . The result is a magnetic field gradient, the sign of 
which has the magnetic force opposing friction everywhere. 

The sign of the velocity gradient in Fig. [12] can be ex- 
plained as follows. The magnetic impulse delivered to a fluid 
particle inside a J front is proportional to SB. The transport 
of field lines away from the fronts reduces SB and hence Svi . 
The result is a velocity gradient with the observed sign. In 
fact there is a feedback mechanism at work: the reduction 
of Svi tends to make the velocity gradient steeper, which 
increases the rate of field line migration (cf. eq. I14p . And 
so on. This explains why jumps in the fluid variables decay 
exponentially rather than (say) linearly with time. 

As noted above, at t — 0.1 magnetic and coUisional 
forces are approximately equal between the advancing front 
and the contact discontinuity, differing by < 1%. This may 
seem surprising, since one might think that balance between 
forces could not occur for t < 1. This is because Tin is ba- 
sically the time it takes for an ion to "lose memory" of its 
initial state through collisions with the neutrals. It is readily 
shown that, in the absence of magnetic forces, the velocity of 
an ion in a uniform one-dimensional flow of neutrals is given 
by Vi{t) = Uioe"*'^^'" +Vn (l — e~*/'"'") , where vio is the ion's 
initial velocity. Approximate balance between magnetic and 
drag forces, however, requires only that the two forces be 
comparable in magnitude, yielding a residual that is much 
smaller than the magnitude of the individual forces them- 
selves. Effective force-balance in the ions behind the J front 
is then possible when the ratio of the net (linearized) accel- 
eration relative to the drag acceleration becomes negligible, 
i.e., when {dvi/dt)[Tin/{vn — «i)] ^ 1 (see eq. [2]). Taking 
dvi,ch/dt ~ Vi^ch/t, where v^^ch is a characteristic post-jump 
ion velocity, and using Un — — in this region (see Fig. 
12), this condition is equivalent to t/rin S> Vi.ch/wn. If the 
magnitude of the ion velocity in the post-shock region is 
sufficiently reduced by the magnetic impulse at the front so 
that ?;i,ch/wn ^ 1, it follows that near force-balance in the 
ions can happen even for i/rin < 1. This is exactly the situa- 
tion that is depicted in Figs. 12 and 14: the magnetic "kick" 
the ions receive at the jump dramatically reduces the ion 
velocity, resulting in Ui,ch/vn ^ 0.05 behind the front, allow- 
ing near-equality of forces in that region at t — 0.1 (> 0.05). 
For this case, the coUisional drag of the streaming neutrals 
on the ions is almost enough to balance the magnetic force 
due to the field gradient behind the front. 

The ion motion is almost but not exactly force-free at 



^ The total flux is just the flux swept up by the J fronts, which 
depends only on Bg and the front velocities. The latter are ztuiAo 
whether or not friction is present. 
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Figure 17. Solid: Ion velocity at t = 1, 3, and 5. Dashed: Neutral 
velocity at t=5. 

t = 0.1. Outside the g-star region, gas drag pulls the ions 
toward the contact discontinuity and the magnetic force 
pushes them away (Fig. I15|l . Because the magnetic force is 
slightly larger in magnitude the ions slow down, to speeds 
which are almost but not quite zero when they pass into 
the g-star region (Fig. I13p . There both forces change sign 
(Fig. I16p : the very small net force still tends to slow the 
ions down, allowing them to come to rest at the contact 
discontinuity. 

Figures [17] and [18] describe the flow at three times of 
order unity, the largest of which is « 0.05 yr. The transport 
of field lines away from the J fronts has attenuated the latter 
significantly, with a corresponding increase in B within the 
magnetic core. These figures show the earliest phase in the 
formation of two multifluid, MHD shock waves by the cloud- 
cloud collision. At t ~ 1, each multifluid shock is comprised 
of a J shock in the neutral fluid with a nascent magnetic 
precursor extending upstream from its J front. Although 
we made some unrealistic assumptions in order to work the 
problem analytically (no radiative cooling, Un independent 
of x), the mismatch between the time here (< 0.1 yr) and 
the time to accelerate the neutral fluid (~ 10* yr) leaves no 
doubt that the solution is qualitatively correct. 

One does not expect the entire (charged-|-neutral) flow 
will approach a steady conflguration until much greater 
times ~ 10* yr. However it is reasonable to suppose, given 
the large "signal speed" in the charged fluid, that the mag- 
netic precursors would reach a quasi-steady state at shorter 
times. This is clearly not the case at t ~ 1 where, for ex- 
ample, the precursors are increasing in width (Fig. I17|l and 
the magnetic flux inside the core is growing (Fig. I18|) . This 
evolving structure is determined mostly by waves propagat- 
ing into the charged fluid. Thus the linear scale (oc t) is set 
by the steady speed of the J fronts. In mathematical terms: 
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Figure 18. Magnetic field perturbation at t = 1, 3, and 5. 

the solution is "mostly" the homogeneous solution, i.e., the 
transient response of the charged fluid to the discontinu- 
ous initial conditions. We refer to this evolutionary phase 
henceforth as the "ion-electron transient." The ion-electron 
transient is unobservable; however it obviously affects sub- 
sequent, possibly observable, phases. 

Figures [19] and [20] displav the flow at three times ~ 10, 
the largest of which corresponds to « 0.5 yr. Now the solu- 
tion is entirely the particular (=driven) solution. The em- 
bedded neutral shocks influence the flow mainly through the 
neutral velocity profile, which dictates the spatial depen- 
dence of the friction force. The fact that the neutral shocks 
are propagating, for example, is relatively unimportant. All 
traces of wave physics have disappeared but the precursors 
are still evolving. The flow in this phase is determined by 
the competition between the advection of field lines inward, 
driven by the ion motions, and the diffusion of field lines 
outward, due to the resulting magnetic field gradient. The 
dominance of diffusion is signalled by the linear scale, which 
now is oc \/^. We refer to this period as the "diffusion phase" 
of multifluid shock formation. 

It would be useful to know what physical effects ter- 
minate the diffusion phase, whether this phase lasts long 
enough to be observable, and whether the charged flow ap- 
proaches a quasi-steady configuration. Unfortunately, the 
solutions for t > 50 cannot be discussed here because non- 
linear effects start to become important (e.g.. Fig. I20|l . Cal- 
culations in progress will use numerical methods to include 
nonlinearity, radiative cooling, and variations in the density 
of the neutral gas. We are also exploring similarity solu- 
tions (for the diffusion phase) of the nonlinear equations of 
motion. If a similarity solution existed, it would rule out 
quasi-steady flow of the charged particles. It might also re- 
veal which physical effects usher in the next evolutionary 
phase. The flow in Fig. [19] and [20] certainly appears to be 
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Figure 19. Solid: Ion velocity at t = 10, 30, and 50. Dashed: 
Neutral velocity at t=50. 




Figure 20. Magnetic field perturbation at t = 10, 30, and 50. 



self-similar. It also satisfies the prerequisites: there are no 
boundary conditions and the system has "forgotten" its ini- 
tial conditions. 



5 SUMMARY 

This paper can be summarized as follows; 

(i) Our objective was to understand basic physics gov- 
erning the formation of multifluid, MHD shock waves from 
plausible initial conditions. We focused on the earliest stages 
of this process, which have not been explored clsowhoro. 

(ii) We treated the plasma as separate fluids of charged 
and neutral particles which are coupled by ion-ricutral fric- 
tion. We exploited the large inertial mismatch between the 
neutral and charged fluids to simplify the calculation. On 
time scales < 10^ yr (tjrpically), the neutral fluid evolves 
as if the charged particles were absent. At sufficiently early 
times the neutral flow can be calculated by solving Euler's 
equations. The flow of the charged fluid is driven by and 
slaved to the neutral flow by friction. 

(iii) We calculated the charged flow for special cases 
where the linearized equations of motion are accurate, and 
carried out an extensive analysis of linear MHD waves driven 
by friction. The physics of driven waves is embodied in cer- 
tain Green functions which describe wave propagation on 
short time scales, ambipolar diffusion on long time scales, 
and transitional behavior at intermediate times. 

(iv) As an illustrative example, we simulated the col- 
lision of two identical clouds with Av = 20 km s~^. The 
simulation incorporated a few unrealistic approximations. 
We have argued that the results are qualitatively correct 
and illustrate the basic physics. Realistic solutions will be 
presented elsewhere. 

(v) We found that the formation of a multifluid shock 
wave proceeds through two initial phases; an "ion-electron 
transient" and a "diffusion phase." In the former, the cloud- 
cloud collision produces J shocks in the neutral ffuid which 
drive wavelike transients into the charged fluid. The tran- 
sients quickly evolve into magnetic precursors on the J 
shocks, wherein the ions undergo force free motion and the 
magnetic field grows steadily in time. In the diffusion phase, 
the charged flow continues to evolve in what appears to be 
self similar fashion. The magnetic precursors do not become 
steady at the largest times wc can study, which are deter- 
mined by the onset of nonlinearity. 
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Figure A2. Contour used to evaluate G+(x,t) for x < t and 
G—{x,t) for X < —t. See text. 



APPENDIX A: FOURIER INTEGRALS FOR 
THE GREEN FUNCTIONS 

In order to find the Green functions it is first necessary to 
calculate G+ {x, t) and G_ (x,t). This involves the evaluation 
of certain Fourier integrals which are almost identical to 
integrals that appear in the Green functions for the Klein- 
Gordon equation. It turns out that "our" integrals can be 
evaluated by exactly the same technique used in the Klein- 
Gordon problem. In i; i]Alf[X2l we closely follow the technique 
of Wyld (1999, see pp. 573fr). 



Al Fourier Integral for G+ 

The function 



47r 



R 



(Al) 



can be evaluated by contour integration as follows. Noting 
that _R — > fc as jfcj oo, we see that 



In the limit \k\ — > oo, the integrand in eq. (|Aip goes to 
zero exponentially in the upper k plane li x — t > Q and in 
the lower k plane ii x — t < 0. It follows that we can find 
G+{x,t) by integrating around the semicircular contour in 
Fig. lAll if X > t or the semicirle in FiE. IA2l if x < t, and then 
taking the limit of infinite radius, p. Because of the branch 
cut in 7? it is important to remember that the horizontal 
part of each semicircle lies "just above" the real k axis (see 
^3:311 . 

Consider first the case x > t (Fig. IA1|) . Since there are 
no singularities in the upper k plane, we find immediately 
that 



G+{x,t)=0 a X > t. 



(A3) 



Now consider the case x < t (Fig. IA2|) . By Cauchy's Theo- 
rem, the contour in Fig. IA2l can be deformed into the one in 
Fig. I A3I without changing the integral. If we simultaneously 
change variables from k to z, where 



k = kc cosh z, 

then expression ([Aip becomes 



(A4) 



G+(x,t) = f e-*/^ 
47r 



dz exp [tkc {x cosh z ~ t sinh z)] , (A5) 



where the contour is now a line segment, as indicated by the 
limits of integration. To proceed it is necessary to distinguish 
whether < i or \x\ > t. 

Suppose first that x < t and < t (i.e., ~t < x < t). 
Introduce the parameter 6 defined implicitly by 

= ■ (A6) 



cosh 9 = — — 
Then 

X — sinh 9, 
t — cosh d, 
where 
ax,t) 



(A7) 
(A8) 

(A9) 



i(kx — Rt) 



^ik{x — t) 



as I A: I 



(A2) 



After eliminating x and t in favor of ^ and 9 in eq. (|A5|I , and 
using some identities for hyperbolic functions, we find that 
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Imz'n 




\x\ — ("cosh-;/) 
and 

t — (^sinhi/), 
where 



(A18) 
(A19) 



Rez' 



({x,t) = ^/x'^-f^. (A20) 

Eliminating {x,t) in favor of ( and ip changes eq. (|XT6)) to 
/.() 

G+(j:,t) = -^e~*/^ / dz expl-ikcC cosh {z + ip)].{A21) 
47r /„ 

27^^ 

Changing variables to z' — z + changes (|A21[) to 



G+(x,t) = ^e-*/^ 



dz' exp (-ifccC cosh 2') (A22) 



Figure A4. Contour used to simplify expression IIAIH (after 
Wyld (1999, see p. 579). 



and considerations analogous to the ones leading from ex- 
pression HA11|) to (|A12|) give 



f dz exp [tk^sinhie - z)]. (AlO) ^+(2:, t) = e y 

J 27ri ^ 



dz' exp (— zfccC cosh 2'). (A23) 



Changing the variable of integration again from 2 to z' 
z — 9 yields 



G+{x,t) = 



« -t/2 

— e ' 
4tt 



dz' exp (— ifcc^ sinh 2') . (All) 



The final transformation 2' = ly yields 
1 

G+{x,t) = — e~*^^ / dy exp {—ikc cosy) . 



(A24) 



2-iTi-e 

The contour for this integral is the left side of the rect- 
angle in Fig. IA4I Now the integral around the whole rectan- 
gle is zero because the integrand has no singularities inside 
the contour. Further, the integrals along the top and bottom 
sides cancel one another because sinh(a; -I- 27ri) = sinhi. It 
follows that 

G+ (a;, t) = - ^ e"*/^ y dz' exp (-ifccC sinh 2') . ( A12) 

If we make one more change of variables such that 2' = 
t{y + 7r/2), we find that 

r.3ir/2 



Noting that the ordinary Bessel function of order zero is 

r.27r 



Jo{x) 



2tt 



dy exp {±ix cos y) . 



we conclude that 



G+{x,t) 



-t/2 



Jo (fccC) if a; < -t. 



G+{x,t) = -. 



-t/2 



dy exp {kc£, COS y) . (A-13) 



To summarize: 

G+{x,t)={ ie-'/2/o(fccO ii-t<x<t 
a x>t 



(A25) 



(A26) 



(A27) 



7r/2 



Now the modified Bessel function of order zero has the in- 
tegral representation 



Jo (a;) = — / dy exp {±xcosy) . 



(A14) 



A2 Fourier Integral for G 

This calculation is very similar to the preceding one; details 
are included for the morbidly curious. We need to evaluate 



Comparing expressions (|A13|I and (|A14|) . and noting the pe- 
riodicity of the integrand in the former, we finally conclude 
that 

G+ (a;, t) = i e"'/' Jo (k^) ii - t < x < t. ( A15) 

Next suppose that x < t and |a:| > t (i.e., x < — t). 
Then x = — |a;| and eq. HA5|) can be written in the form 



G_(.,t) = ^e-/^ 



^fcg.(fcx + flt) 

R 



Now the exponential factor behaves like 



as Ifcl 



(A28) 



(A29) 



G+{x,t) = - — e~*^^ / d2 exp [—ifcc (|a;| cosh 2 + t sinh 2)] 
47r /„ 

J 2m 



so we use the contour in Fig. lAll if x > and the one in 
Fig. IA2I if a; < —t. Since there are no singularities in the 
upper k plane we find immediately that 



(A16) 



(A30) 



Now define tp by 
cosh-i/) ; ' ' 



Vx^ - t2 ' 



(A17) 



so that 



G-{x,t)^0 ifx>-t. 

U X < —t we change variables from fc to 2 (cf. eq IA4|I and 
deform the contour to the one in Fig. IA3I with the result 

G-(a;,t) = — e"'''^ / d2 exp [-ifccC cosh (2 - e)] .(A31) 

J2m 
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Comparing the integrals in eq. HA16|I and (|A31|I , and noting 
that the former gives a result which is independent of 0, we 
see that 



{x, t) = ^ e-'l^ Jo (fccC) ifa; < -t. 
To summarize: 



(A32) 



i e~*/2 Jo (fc,C) if s < 
G-(x,t)={ if-t<x<t . (A33) 

a x>t 



APPENDIX B: ASYMPTOTIC BEHAVIOR OF 
THE GREEN FUNCTIONS 

For times < 1 the Green functions exhibit complex behavior 
including aspects of both wave propagation and diffusion 
(Fig. |4H6]). However the properties of the dispersion rela- 
tion f ij3.3p suggest that diffusion dominates at large times, 
and the Green functions should behave accordingly. We now 
show that this is indeed the case by calculating the Green 
functions for t ^ 1. 

First consider the integral for G+{x,t) in eq. (|A1[I . The 
time dependent part of the integrand is 



-t/2 -iRt 

■ e ' e 



(Bl) 



and this behaves differently for large and small wave num- 
bers. Since R{k) is real for > fee, we can neglect contri- 
butions to the integral from k > kc = 1/2. For < 1/2 the 
definition of R [eq. (|34|) ] says 



i?(fc) = ^vr^4p (ifei<i/2) 



(B2) 



Since this factor decays faster than exp (—t/2) for all wave 
numbers, we find 



G-{x,t)^0 ift>l 

and taking the difference of G+ and G- gives 
1 



G{x,t) 



/47rf 



: e « ff t > 1. 



(BIO) 



(Bll) 



Taking partial derivatives gives the other Green functions: 



G{x,t) X A{x,t) : 



Y x ^ 




\2t, 


' 2t 



G{x,t) if t > 1, (B12) 



;'(x,t)«r(x,t)^-(j) G(a;,t) ift 



and 
G' 

Notice that 

\G{x,t)\ < \G{x,t)\ 
and 

\G\x,t)\ < t-^'^ \G{x,t)\ 
for large times. 



> 1. 



(B13) 
(B14) 
(B15) 



6"*"+' =exp|-| 1- (l-4fc^)'''' I (|fc| < 1/2). (B3) 

As expected, the integral is dominated by contributions from 
k <^ 1 (i.e., long wavelengths) at large times. Setting 



e + Ri e 
and 

inside the integral sign, we find 



G+{x,t)^- 



dk exp (^—k^t + ikx^ (t 3> 1). 

The integral is easily evaluated to give 
G+{x,t) ^ —^e~'^ if t > 1. 



1 



(B4) 
(B5) 

(B6) 
(B7) 



Next consider the integral in expression (|A28|I for 
G-{x,t). Now the time dependent factor in the integrand 
is 



-t/2 +iRt 



(B8) 

Once again contributions to the integral from |fc| > 1/2 can 
be neglected. For |fc| < 1/2, the definition of R(k) implies 



= exp ■ 



{-f [l+(l-4fc')'^']} (|fe|<l/2).(B9) 



